#!/usr/bin/python
# plot_logN.py

import math
import numpy
from matplotlib import use as useBackend
useBackend('Agg')
import matplotlib.pyplot as ppt


def cheat111():
    '''Cheat to clear the plot.'''
    ppt.subplot(311)
    ppt.plot([0,1],[1,2])
    ppt.subplot(131)
    ppt.plot([0,1],[1,2])
    ppt.subplot(111)
    ppt.hold(True)

colorlist = ['k', 'b', 'g', 'r', 'c', 'm', 'y', 'indigo', 'darkorange', 'saddlebrown', 'forestgreen', 'deeppink', 'slategray', 'maroon', 'lime', 'bisque', 'aquamarine', 'darkcyan', 'darkmagenta', 'mediumpurple', 'olive', 'purple', 'salmon'] #23

def choosecolor(j, colorlist):
    ncolors = len(colorlist)
    while True:
        if j > (ncolors - 1):
            j = j - ncolors
        else:
            colornow = colorlist[j]
            break
    return colornow

plottinglimits = numpy.array([[11.5,12.5],[13.5,14.5],[15.5,16.5],[17.5,18.5],[19.5,20.5]])


def plot(simpars, phypars, compN, N, allskyN):
    '''Creates a plot similar to Han's figure 1.

    REQUIRES inputs constructed using the 'vanilla' limitset.'''
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    # Relevant distances, conversions; note N, N_bM are already in log deg**-2
    ba = blocs*PtScaleba*180./math.pi
    plottingba = numpy.zeros(2*bres)
    plottingba[bres:] = ba[:]
    plottingba[:bres] = -ba[::-1]
    # Plot for Mv = 12,14,16,18,20
    cheat111()
    legendlist = []
    plotlist = []
    index = 0
    N_bM_index = 0
    for j in range(0, len(plottinglimits)):
        hiplotlim = plottinglimits[j,1]
        loplotlim = plottinglimits[j,0]
        thisarray = numpy.zeros(bres)
        thisfinalarray = numpy.zeros(2*bres)
        for i in range(0, len(limitset)):
            lolim = limitset[i,0]
            hilim = limitset[i,1]
            limdiff = hilim - lolim
            # If it falls entirely within the range
            if (hilim <= hiplotlim) and (lolim >= loplotlim):
#                print "Case a", hiplotlim, loplotlim, hilim, lolim
                thisaddition = 10**N[i,:]
                thisarray = thisarray + thisaddition
            # If it overlaps the high edge of the range
            if (hilim >= hiplotlim) and (lolim > loplotlim) and (lolim < hiplotlim):
#                print "Case b", hiplotlim, loplotlim, hilim, lolim
                overlapwidth = float(hiplotlim - lolim)
                fullwidth = float(hilim - lolim)
                thisaddition = (10**N[i,:])*overlapwidth/fullwidth
                thisarray = thisarray + thisaddition
            # If it overlaps the low edge of the range
            if (hilim < hiplotlim) and (lolim <= loplotlim) and (hilim > loplotlim):
#                print "Case c", hiplotlim, loplotlim, hilim, lolim
                overlapwidth = float(hilim - loplotlim)
                fullwidth = float(hilim - lolim)
                thisaddition = (10**N[i,:])*overlapwidth/fullwidth
                thisarray = thisarray + thisaddition
            # If it entirely encompasses the range
            if (hilim > hiplotlim) and (lolim < loplotlim):
#                print "Case d", hiplotlim, loplotlim, hilim, lolim
                overlapwidth = float(hiplotlim - loplotlim)
                fullwidth = float(hilim - lolim)
                thisaddition = 10**N[i,:]*overlapwidth/fullwidth
                thisarray = thisarray + thisaddition
        thisfinalarray[bres:] = thisarray[:]
        thisfinalarray[:bres] = thisarray[::-1]
#        print thisfinalarray
        # And plot this element
        legendname = str(loplotlim)+' < Mv < '+str(hiplotlim)
        legendlist.append(legendname)
        cn = choosecolor(index, colorlist)
        thisplot = ppt.semilogy(plottingba, thisfinalarray, c=cn)
        plotlist.append(thisplot)
#        ppt.semilogy(-ba, 10**N[i,:], c=cn)
        ppt.semilogy(bvals, 10**N_bM[N_bM_index,:], c=cn, ls='none', marker='o')
        ppt.semilogy(-bvals, 10**N_bM[N_bM_index,:], c=cn, ls='none', marker='o')
        index = index + 1
        N_bM_index = N_bM_index + 2
    ppt.axis([-91, 91, 10, 10**5])
    ppt.xlabel("Galactic latitude")
    ppt.ylabel("Averaged stars deg$^{-2}$ dM$_{V}$$^{-1}$")
    ppt.title("Measured vs. simulated stellar number densities")
    ppt.legend(plotlist, legendlist)
    rhostr = str(rho0)[:6]
    n = len(rhostr)
    n = 6 - n
    rhostr = rhostr + "0"*n
    cstr = str(c0)[:7]
    n = len(cstr)
    n = 7 - n
    cstr = cstr + "0"*n
    intstr = str(1./float(intervals))[:5]
    n = len(intstr)
    n = 5 - n
    intstr = intstr + "0"*n
    plotname = outdir+"Search/HanLikeFig1_rho_"+rhostr+"_c_"+cstr+"_int_"+intstr+".pdf"
    ppt.savefig(plotname)
